(fig1i_vt <- read_tsv(here("data-raw/tables/Fig1i_venntable.tsv"))) %>% 
  kbl() %>%
  kable_material(c("striped", "hover"))

(dna_lines_vt <- read_tsv(here("data-raw/tables/DNALines_venntable.tsv")))%>% 
  kbl() %>%
  kable_material(c("striped", "hover"))

(dna_lines_pertarg_vt <- read_tsv(here("data-raw/tables/pertarget_venntable.tsv")))%>% 
  kbl() %>%
  kable_material(c("striped", "hover"))

(nnn04_nnn05_vt <- read_tsv(here("data-raw/tables/nnn04m_nnn05m_venn_table_df.tsv")))%>% 
  kbl() %>%
  kable_material(c("striped", "hover"))

## enter labs here to exclude from study -- will make it easier to switch between common chunks of code
lab_disregard <- c("NNN02m", "NNN03m", "NNN04m", "SSN01m", "SSN02m")

## count the number of total genommic coordinates for each target site for each lab
total_genomic_positions_1 <- off_target_df %>% 
  filter(!(lab %in% lab_disregard)) %>%
  group_by(target_site, lab) %>% 
  ## including the genomic coordinate col gave me a "Genomic Coordinate" col so excluding 
  count(target_site, lab,  name = "total_number_of_sites") %>%
  ## verified numbers by counting rows of identified_matched.txt files for each target site/lab
  arrange(target_site, lab) %>%
  left_join(lab_sample_df)

ggplotly(ggplot(total_genomic_positions_1) + 
    geom_bar(aes(x = target_site, y = total_number_of_sites, fill = lab), 
             position = position_dodge(preserve = "single"), stat = "identity") +
    scale_y_continuous(labels = comma) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Total Number of Sites", 
         fill = "Lab") + scale_fill_manual(values=cbPalette))

total_genomic_positions_1 %>%
  kbl() %>%
  kable_material(c("striped", "hover"))

lab_disregard <- c("NNN02m", "NNN03m", "NNN04m", "SSN01m", "SSN02m", "NSN02m", "NNN05m")

## identify sites that are found in both replicates 
#################### DF CURATION ############################
## derivate of previous reports' variables' "counts_by_lab"
nnn01_nsn01_overlap_sites <- summary_df %>%
  filter(!(lab %in% lab_disregard)) %>%
  # Only including off target edits
  filter(!Site_Substitution_Number == 0) %>% 
  ## Excluding redundant and unnecessary columns
  # Extra columns results in non-matching counts by lab
  select("Genomic Coordinate", target_site, lab, Nuclease_Read_Count) %>% 
  pivot_wider(names_from = "lab",
              values_from = "Nuclease_Read_Count", 
              values_fill = 0) %>%
  ## NNN0 vs. NNN1: Concordant Off-Target Coordinates
  filter(NNN01m != 0, NSN01m != 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(NNN01m, NSN01m)),
         SD = sd(c(NNN01m, NSN01m)),
         CV = SD / AVG) %>%
  arrange(desc(CV))

#nnn01_nsn01_overlap_sites %>%
#  kbl() %>%
#  kable_material(c("striped", "hover"))

nnn01_nsn01_overlap_sites_plot <- nnn01_nsn01_overlap_sites %>% 
  rename(genomic_coordinate = "Genomic Coordinate") %>%
  select(genomic_coordinate, NNN01m, NSN01m, target_site) 
#  pivot_longer(cols = c("NNN01m", "NSN01m"), names_to = "lab", values_to = "read_count") %>%
#  ungroup() %>%
#  arrange("read_count")

#ggplotly(ggplot(data=nnn01_nsn01_overlap_sites_plot,aes(x=read_count,fill=lab)) + 
#  geom_bar(data=subset(nnn01_nsn01_overlap_sites_plot,lab=="NNN01m")) + 
#  geom_bar(data=subset(nnn01_nsn01_overlap_sites_plot,lab=="NSN01m"),aes(y=..count..*(-1))) + 
#  scale_y_continuous(breaks=seq(-40,40,10),labels=abs(seq(-40,40,10))) + 
#  coord_flip())

ggplotly(ggplot(nnn01_nsn01_overlap_sites_plot) + 
   geom_point(aes(x = NNN01m, y = NSN01m), alpha = 0.25) + 
     scale_x_log10() + 
  scale_y_log10() + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = "NNN01m Read Counts", 
       y = "NSN01m Read Counts"))

lab_disregard <- c("NNN02m", "NNN03m", "NNN04m", "SSN01m", "SSN02m", "NNN01m", "NSN01m")

nnn05_nsn02_overlap_sites <- summary_df %>%
  filter(!(lab %in% lab_disregard)) %>%
  # Only including off target edits
  filter(!Site_Substitution_Number == 0) %>% 
  ## Excluding redundant and unnecessary columns
  # Extra columns results in non-matching counts by lab
  select("Genomic Coordinate", target_site, lab, Nuclease_Read_Count) %>% 
  pivot_wider(names_from = "lab",
              values_from = "Nuclease_Read_Count", 
              values_fill = 0) %>%
  ## Concordant Off-Target Coordinates
  filter(NNN05m != 0, NSN02m != 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(NNN05m, NSN02m)),
         SD = sd(c(NNN05m, NSN02m)),
         CV = SD / AVG) %>%
  arrange(desc(CV))

#nnn05_nsn02_overlap_sites %>%
#  kbl() %>%
#  kable_material(c("striped", "hover"))

nnn05_nsn02_overlap_sites_plot <- nnn05_nsn02_overlap_sites %>% 
  rename(genomic_coordinate = "Genomic Coordinate") %>%
  select(genomic_coordinate, NNN05m, NSN02m, target_site) 

ggplotly(ggplot(nnn05_nsn02_overlap_sites_plot) + 
   geom_point(aes(x = NNN05m, y = NSN02m), alpha = 0.25) + 
     scale_x_log10() + 
  scale_y_log10() + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = "NNN05m Read Counts", 
       y = "NSN02m Read Counts"))

#(nnn05_nsn02_overlap_sites)

lab_disregard <- c("NNN02m", "NNN01m", "NNN05m", "NSN01m", "NSN02m", "SSN01m", "SSN02m")

## count the number of total genommic coordinates for each target site for each lab
total_genomic_positions_2 <- off_target_df %>% 
  filter(!(lab %in% lab_disregard)) %>%
  group_by(target_site, lab) %>% 
  ## including the genomic coordinate col gave me a "Genomic Coordinate" col so excluding 
  count(target_site, lab,  name = "total_number_of_sites") %>%
  ## verified numbers by counting rows of identified_matched.txt files for each target site/lab
  arrange(target_site, lab)  %>%
  left_join(lab_sample_df)

ggplotly(ggplot(total_genomic_positions_2) + 
    geom_bar(aes(x = target_site, y = total_number_of_sites, fill = lab), 
             position = position_dodge(preserve = "single"), stat = "identity") +
    scale_y_continuous(labels = comma) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Total Number of Sites", 
         fill = "Lab") + scale_fill_manual(values=cbPalette))

total_genomic_positions_2 %>%
  kbl() %>%
  kable_material(c("striped", "hover"))

lab_disregard <- c("NNN02m", "NNN01m", "NNN05m", "NSN01m", "NSN02m", "SSN01m", "SSN02m")

## identify sites that are found in both replicates 
#################### DF CURATION ############################
## derivate of previous reports' variables' "counts_by_lab"
nnn03_nnn04_overlap_sites <- summary_df %>%
  filter(!(lab %in% lab_disregard)) %>%
  # Only including off target edits
  filter(!Site_Substitution_Number == 0) %>% 
  ## Excluding redundant and unnecessary columns
  # Extra columns results in non-matching counts by lab
  select("Genomic Coordinate", target_site, lab, Nuclease_Read_Count) %>% 
  pivot_wider(names_from = "lab",
              values_from = "Nuclease_Read_Count", 
              values_fill = 0) %>%
  ## NNN0 vs. NNN1: Concordant Off-Target Coordinates
  filter(NNN03m != 0, NNN04m != 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(NNN03m, NNN04m)),
         SD = sd(c(NNN03m, NNN04m)),
         CV = SD / AVG) %>%
  arrange(desc(CV))

#nnn03_nnn04_overlap_sites %>%
#  kbl() %>%
#  kable_material(c("striped", "hover"))

nnn03_nnn04_overlap_sites_plot <- nnn03_nnn04_overlap_sites %>% 
  rename(genomic_coordinate = "Genomic Coordinate") %>%
  select(genomic_coordinate, NNN03m, NNN04m, target_site) 

ggplotly(ggplot(nnn03_nnn04_overlap_sites_plot) + 
   geom_point(aes(x = NNN03m, y = NNN04m), alpha = 0.25) + 
     scale_x_log10() + 
  scale_y_log10() + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = "NNN03m Read Counts", 
       y = "NNN04m Read Counts"))

#(nnn03_nnn04_overlap_sites)

lab_disregard <- c("NNN02m", "NNN01m", "NNN03m", "NSN01m", "NSN02m", "SSN01m", "SSN02m")

## count the number of total genommic coordinates for each target site for each lab
total_genomic_positions_3 <- off_target_df %>% 
  filter(!(lab %in% lab_disregard)) %>%
  group_by(target_site, lab) %>% 
  ## including the genomic coordinate col gave me a "Genomic Coordinate" col so excluding 
  count(target_site, lab,  name = "total_number_of_sites") %>%
  ## verified numbers by counting rows of identified_matched.txt files for each target site/lab
  arrange(target_site, lab)  %>%
  left_join(lab_sample_df)

ggplotly(ggplot(total_genomic_positions_3) + 
    geom_bar(aes(x = target_site, y = total_number_of_sites, fill = lab), 
             position = position_dodge(preserve = "single"), stat = "identity") +
    scale_y_continuous(labels = comma) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Total Number of Sites", 
         fill = "Lab") + scale_fill_manual(values=cbPalette))

total_genomic_positions_3 %>%
  kbl() %>%
  kable_material(c("striped", "hover"))

lab_disregard <- c("NNN02m", "NNN01m", "NNN03m", "NSN01m", "NSN02m", "SSN01m", "SSN02m")

## identify sites that are found in both replicates 
#################### DF CURATION ############################
## derivate of previous reports' variables' "counts_by_lab"
nnn05_nnn04_overlap_sites <- summary_df %>%
  filter(!(lab %in% lab_disregard)) %>%
  # Only including off target edits
  filter(!Site_Substitution_Number == 0) %>% 
  ## Excluding redundant and unnecessary columns
  # Extra columns results in non-matching counts by lab
  select("Genomic Coordinate", target_site, lab, Nuclease_Read_Count) %>% 
  pivot_wider(names_from = "lab",
              values_from = "Nuclease_Read_Count", 
              values_fill = 0) %>%
  ## NNN0 vs. NNN1: Concordant Off-Target Coordinates
  filter(NNN05m != 0, NNN04m != 0) %>%
  rowwise() %>%
  mutate(AVG = mean(c(NNN05m, NNN04m)),
         SD = sd(c(NNN05m, NNN04m)),
         CV = SD / AVG) %>%
  arrange(desc(CV))

#nnn05_nnn04_overlap_sites %>%
#  kbl() %>%
#  kable_material(c("striped", "hover"))

nnn05_nnn04_overlap_sites_plot <- nnn05_nnn04_overlap_sites %>% 
  rename(genomic_coordinate = "Genomic Coordinate") %>%
  select(genomic_coordinate, NNN04m, NNN05m, target_site) 

ggplotly(ggplot(nnn05_nnn04_overlap_sites_plot) + 
   geom_point(aes(x = NNN04m, y = NNN05m), alpha = 0.25) + 
     scale_x_log10() + 
  scale_y_log10() + 
  facet_wrap(~target_site, ncol = 1) + 
  theme_bw() + 
  labs(x = "NNN04m Read Counts", 
       y = "NNN05m Read Counts"))

#(nnn05_nnn04_overlap_sites)

## str_c joins character vectors via collapse flag. 
## by unique lab sorting we are joining labs into unique combinations
## every represented lab set for each target site and genomic coordinate is present
lab_disregard <- c("NNN02m", "NNN01m", "NNN03m", "NNN04m", "NNN05m")
off_target_lab_set <- off_target_df %>%
  filter(!(lab %in% lab_disregard)) %>%
#  filter(target_site == "E_CTLA4_site_9") %>%
  select(`Genomic Coordinate`, target_site, lab, Nuclease_Read_Count) %>%
  group_by(`Genomic Coordinate`, target_site) %>% 
  summarise(lab_set = str_c(sort(unique(lab)), collapse = "-"))

venn_df <- off_target_df %>%
  select(`Genomic Coordinate`, target_site, lab, Nuclease_Read_Count) %>%
  filter(!(lab %in% lab_disregard)) %>%
#  filter(target_site == "E_CTLA4_site_9") %>%
  pivot_wider(names_from = "lab", values_from = "Nuclease_Read_Count", values_fill = NA) %>%
  left_join(off_target_lab_set)


ggplotly(ggplot(data = venn_df) + 
  geom_bar(aes(fct_rev(as_factor(target_site)), fill = lab_set), 
           # position = "fill"
           position = "fill"
           ) + 
  coord_flip() + 
  theme_bw() + 
  labs(x = "Target Site", 
       y = "Proportion of Sites", 
       fill = "Lab Combinations") + scale_fill_manual(values=cbPalette))

####### tables

## count sites in lab-sets
(num_sites <- off_target_lab_set %>%
  group_by(target_site) %>%
  count(lab_set)) %>% 
  rename("total_off-target_sites" = n) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))

## for number of sites each bin = total # times particular combination shows up 


## read counts 
(read_range_df <- venn_df %>%
  pivot_longer(cols = c("NSN01m", "NSN02m", "SSN01m", "SSN02m"), 
               names_to = "sample", values_to = "read_counts", values_drop_na = TRUE) %>%
    group_by(target_site, lab_set) %>%
  summarise(min_read_count = min(read_counts), 
            max_read_count = max(read_counts))) %>% 
  rename(target_site = target_site)%>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>% 
  kable_styling(font_size = 8)

## curate the lower sequencing depth df
low_seq_disc <- nnn04_nnn05_discordants %>%
  pivot_longer(cols = c("NNN04m", "NNN05m"), names_to = "lab", values_to = "low_seq_discordant") %>%
  filter(!is.na(low_seq_discordant)) %>%
  select(Genomic_Coordinate, target_site, low_seq_discordant)
## low_seq_discordant = read count value of either NNN04m or NNN05m 

## combine the low_seq_disc sites with the sites for the manuscript samples
low_seq_disc_both_reps_df <- left_join(low_seq_disc, dna_lines_coords, by = "Genomic_Coordinate") %>%
  rename(high_seq_rep1 = "SSN01m", high_seq_rep2 = "SSN02m")

## create venn df's 
low_seq_lab_set <- low_seq_disc_both_reps_df %>%
  select(Genomic_Coordinate, target_site.x, low_seq_discordant, high_seq_rep1, high_seq_rep2) %>%
#  mutate_all(~replace(., is.na(.), 0)) %>%
  pivot_longer(cols = c("low_seq_discordant", "high_seq_rep1", "high_seq_rep2"), 
               names_to = "lab", values_to = "Nuclease_Read_Count", values_drop_na = TRUE) %>%
  select(Genomic_Coordinate, lab, Nuclease_Read_Count, target_site.x) %>%
  group_by(Genomic_Coordinate, target_site.x) %>% 
  summarise(lab_set = str_c(sort(unique(lab)), collapse = "-"))

low_seq_venn_df <- low_seq_disc_both_reps_df %>%
  select(Genomic_Coordinate, low_seq_discordant, high_seq_rep1, high_seq_rep2, target_site.x) %>%
  left_join(low_seq_lab_set, by = "Genomic_Coordinate")

## plot
ggplotly(ggplot(data = low_seq_venn_df) + 
  geom_bar(aes(fct_rev(as_factor(target_site.x.x)), fill = lab_set), 
           position = "fill"
           ) + 
  coord_flip() + 
  theme_bw() + 
  labs(x = "Target Site", 
       y = "Proportion of Sites", 
       fill = "Sample IDs") + scale_fill_manual(values=cbPalette))

## for number of sites each bin = total # times particular combination shows up 

## count sites in lab-sets
(low_high_num_sites <- low_seq_lab_set %>%
  group_by(target_site.x) %>%
  count(lab_set)) %>% 
  rename("total_off-target_sites" = n) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))

## for number of sites each bin = total # times particular combination shows up 


## read counts 
(low_high_read_range_df <- low_seq_venn_df %>%
  pivot_longer(cols = c("low_seq_discordant", "high_seq_rep1", "high_seq_rep2"), 
               names_to = "sample", values_to = "read_counts", values_drop_na = TRUE) %>%
    group_by(target_site.x.x, lab_set) %>%
  summarise(min_read_count = min(read_counts), 
            max_read_count = max(read_counts))) %>% 
  rename(target_site = target_site.x.x)%>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>% 
  kable_styling(font_size = 8)